home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Shareware Grab Bag
/
Shareware Grab Bag.iso
/
090
/
byt1186b.arc
/
RUNGKUT.LBR
/
PREDICT.FOR
< prev
next >
Wrap
Text File
|
1986-04-11
|
2KB
|
76 lines
$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c***********************************************************************
subroutine rkp(x,ys,h,neqn,imeth,kt,est,yold,ynew,w,func)
implicit double precision (a-h,k,p-z)
external func
dimension w(13,neqn),kt(neqn),ys(neqn),est(neqn)
dimension al(12),b(12,12),yold(neqn),ynew(neqn)
common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
& a12,a13,er1,er3,er4,er5,er6,er7,er8,
& er9,er10,er11,er12,er13,inum
xs=x
do 5 i=1,neqn
ys(i)=yold(i)
5 continue
call func(xs,ys,kt,neqn)
do 7 i=1,neqn
w(1,i)=kt(i)
7 continue
do 20 j=2,inum
j1=j-1
xs=x+h*al(j1)
do 15 i=1,neqn
ksum=0.d0
do 10 m=1,j1
ksum=ksum+b(m,j1)*w(m,i)
10 continue
ys(i)=yold(i)+h*ksum
15 continue
call func(xs,ys,kt,neqn)
do 17 i=1,neqn
w(j,i)=kt(i)
17 continue
20 continue
c**** evaluate YNEW[i] and the error estimates EST[i]
if(imeth.EQ.1) then
do 30 i=1,neqn
ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
& +a5*w(5,i)+a6*w(6,i))
est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
& +er6*w(6,i))
30 continue
end if
if(imeth.EQ.2) then
do 40 i=1,neqn
ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
& +a5*w(5,i)+a7*w(7,i)
& +a8*w(8,i))
est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
& +er6*w(6,i)+er7*w(7,i)+er8*w(8,i))
40 continue
end if
if(imeth.EQ.3) then
do 50 i=1,neqn
ynew(i)=yold(i)+h*(a1*w(1,i)+a6*w(6,i)+a7*w(7,i)
& +a8*w(8,i)+a9*w(9,i)
& +a12*w(12,i)+a13*w(13,i))
est(i)=h*(er1*w(1,i)
& +er6*w(6,i)+er7*w(7,i)+er8*w(8,i)+er9*w(9,i)
& +er10*w(10,i)+er11*w(11,i)+er12*w(12,i)
& +er13*w(13,i))
50 continue
end if
return
end